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Summary 


The method based on Fourier functional analysis and indicial 
formulation for aerodynamic modeling as proposed by Chin and Lan is 
extensively examined and improved for the purpose of general applications 
to realistic airplane configurations. Improvement is made to automate the 
calculation of model coefficients, and to evaluate more accurately the 
indicial integral. Test data of large angle-of-attack ranges for two different 
models, a 70-deg. delta wing and an F-18 model, are used to further verify 
the applicability of Fourier functional analysis and validate the indicial 
formulation. The results show that the general expression for harmonic 
motions throughout a range of k is capable of accurately modeling the 
nonlinear responses with large phase lag except in the region where an 
inconsistent hysteresis behavior from one frequency to the other occurs. The 
results by the indicial formulation indicate that more accurate results can 
be obtained when the motion starts from a low angle of attack where 
hysteresis effect is not important. 
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1. INTRODUCTION 


Due to the rapid development of on-board computers, a new air 
combat tactic, called supermaneuverability, has become more feasible and 
realistic than before. Fighter designers have recognized the applications of 
supermaneuverability as one of the major features of next-generation 
fighters. With the advantage of simplicity, flying up to and beyond post- 
stall(PST) region has been considered as one way of achieving 
supermaneuverability(Ref.l). One of the applicable PST maneuverings is a 
180-degree change of heading with the additional constraint of returning 
to the point of departure at the initial speed and altitude (Fig.l). However, 
the usefulness of PST maneuvering in terms of tactical purpose still is in 
question. Instead of conducting flight tests, such as X-31 demonstrator, 
flight simulation is an easier and more flexible way to evaluate the 
advantage of PST maneuverings. 

In flight simulation, aerodynamic forces and moments are needed. 
According to Tobak and Schiff in Ref.2, the major difficulty in calculating 
the nonlinear aerodynamic forces and moments acting on a rapidly 
maneuvering aircraft is that these airloads are, in general, determined not 
only by the instantaneous motion variables, for example a and a , but also 
by all of the prior states of the motion up to the current state. As a result, 
the currently-used linear and locally-linearized quasi-steady aerodynamic 
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models incorporated in the look-up tables of wind tunnel test data have 
become increasingly inadequate for the simulation of a rapid maneuvering 
flight. Although a simple empirical surface pressure model was proposed 
lately in Ref 3. to improve the accuracy of nonlinear aerodynamic modeling 
for a 65-deg. delta wing subjected to large-amplitude high-rate oscillation 
in roll, it is not applicable to general configurations. Therefore, a more 
general aerodynamic modeling technique is needed. 

In applications of linear potential flow theory, aeroelasticity 
researchers(Ref.4 and 5) utilized Fourier Transform to relate the 
aerodynamic responses of a step change in the angle of attack of a wing to 
those of harmonic oscillatory motions. The transient aerodynamic response 
corresponding to a step change in the angle of attack, called an "indicial 
function", has been calculated for several classes of isolated wings(Ref.4-7). 
Further applications of the indicial function were considered by Tobak and 
his coworkers by extending the concept of indicial function into the 
nonlinear aerodynamic regimes(Ref.2 and 8). In addition, the method of 
separating the time-history data into in-phase and out-phase components 
has been successfully carried out for the type of response with small- 
amplitude oscillations(Ref.9). However, the above-mentioned simple model 
which only includes the fundamental frequency and small-amplitude is not 
applicable to nonlinear aerodynamic responses involving dynamic stall and 
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vortex lag. 

For this reason, Chin and Lan(Ref.lO) proposed a general model 
based on Fourier analysis to analyze the force and moment data obtained 
in large amplitude forced oscillation tests at high angles of attack. Test 
data for a 70-deg. delta wing were used to verify this method of analytically 
modeling responses of harmonic motions at different reduced frequencies. 
In addition, harmonic ramp motions for the 70-deg. delta wing also were 
calculated to verify the indicial formulation in their paper. 

Since the method of Fourier functional analysis uses a two-level 
gradient method to determine the model coefficients for all linear and 
nonlinear terms in the general model, questions arise as to the uniqueness 
of results. In other words 

•are the model coefficients sensitive to the initial guesses? 

•If the solutions are sensitive to initial guesses, are the results 
obtained with different initial guesses acceptable? 

•If not all results by different initial guesses are acceptable, what 
would the constraint and criteria be to determine the "best" set of 
coefficients? 

In this research, two more sets of test data obtained at NASA 
Langley Research Center will be used to conduct a more extensive testing 
on the method of Fourier functional analysis as well as the indicial 
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formulation. The results have shown that an outer optimization procedure 
with constraints should be added to the original Fourier functional analysis 
due to the complexity of the present highly-nonlinear optimization problem 
with constraints. In addition, the results in validating the indicial 
formulation also show that some special treatments should be made in the 
practical applications of the indicial formulation. 
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2. FOURIER FUNCTIONAL ANALYSIS 

2.1 Theoretical Model 

In the Chin and Lan’s study, a general formulation for a response C L 
at time t can be written as 


C L (t ) =C L (0) + zero-lag response 

+ Tcl [t-x; a(x ) , Cc(x) ] da J —dx 
Jo a dx 

+ — pCr [t-x / oc(x 
vjo L * 


(i: 


&(x ) ] d( * (T idx 
dx 


Here, the zero-lag response represents the virtual mass effect in 2-D 
incompressible flow which includes the effect of a and is identical in every 
harmonic motion. Therefore, it is independent of the time history of motion. 
The last two terms involving the time integrals represent the summation 
of indicial responses at time t due to changes in a and a in the prior 
motions. The key task to determine the time integrals in Eq.(l) is to find 
an analytical form for C L in terms of a(t) and a (t). Then, the time response 
at a given time t can be calculated through the integrals by substituting the 
derivatives of C L . 

According to the linear theory of unsteady aerodynamics, an 
unsteady aerodynamic response can be separated into a product of an 
amplitude function and a phase function in harmonic motion. The 
amplitude function is a function of motion variables as well as their time 
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rate of change. On the other hand, the phase function is a function of 
frequency and depends on the phase lag between the response and the 
excitation. For 2-D incompressible flow, Theodorsen introduced a phase 
function(Ref.4 and 5) which can be determined numerically by analyzing 
the responses obtained at different frequencies with the same amplitude in 
harmonic oscillations. Chao and Lan(Ref.ll) used this approach successfully 
to calculate the indicial lift function for a plunging wing. 

Based on a similar idea, a more general model which is applicable to 
nonlinear aerodynamic responses involving dynamic stall and vortex lag is 
formulated to be a sum of the products of amplitude functions and phase 
functions for frequencies which are multiples of the test frequency in 
harmonic motion, instead of only taking the response to the test frequency 
as in the linear theory. That is, 

C L =C 0 +E (amplitude function) j * (phase function) j (2) 

j 

For a harmonic motion at a given frequency, the response can be 
decomposed in terms of k and t’ by Fourier-analyzing this response over one 
period. The complete expression is written as 


C l =Ao+Ai c °s ( kt 7 ) +A 2 cos ( 2 kt ; ) +A3COS ( 3 kt / ) 
+ 133310 (kt^ ) +B 2 sin ( 2 kt / ) +B 3 sin ( 3 kt / ) 
+ ; . . 


(3a) 
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(3b) 


a l = oc m + a o cos ( kt/ ) 

a=a D cos(kt / ) (3c) 

Cc= (-a 0 k) sin (kt/ ) (3d) 

Here, k is the reduced frequency, t’ is the nondimensional time, is the 
mean angle of attack, and a 0 is the amplitude of angle-of-attack change. 
Following the convention in the classical airfoil theory, the analysis is best 
performed in complex form as follows: 

C l =A 0+ (Ai-iBi) e ikt ' + (A 2 -iB 2 ) e i2kt/ (4a) 

+ (A 3 -iB 3 ) e i3kt/ + . . . 

a=a 0 e ikt ’ (4b) 

a =(ia Q k) e lkt (4c) 

An important step herein is to convert Eq.(4) into a formula in terms of a(t’) 
and a (t’). A "successive Fourier analysis" has been used effectively to split 
cos(n0) and sin(n0) into n+1 terms which include cos n 0 and sin n 0 and other 
cross-product terms. Once the coefficients Aj,Bj at different frequencies are 
obtained by Fourier analysis, the next question is what the appropriate 
expression is for the phase function to represent the lag effects throughout 
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the whole range of k. From past experience(Ref.lO and 12), it was found 
that Pade approximants provide an accurate approximation of the 
theoretical phase functions. Therefore, in the present model the response 
will be put in the following form including the products of amplitude 
functions and phase functions: 

— ^0 (k) 

+ c x * [E n 0c + E 21 d + (H n a + H 21 a ) * (1 - PD X ) ] 

+ c 2 * [Ei 2 0t 2 + e 22 0c 2 + (H 12 a 2 + H 22 aa + h 32 & 2 ) (5) 

* (1 - PD 2 ) ] 

+ C 3 ♦ + E 23 d 3 + (H^^OC 2 + H 2 3 (X Cl + H 33 OCOC + E^ 3 d ) 

* (1 - PD 3 ) ] 

+ . . . 

where PD’s are Pade approximants of order 2 and are defined as 

Pi-j (ik) 2 + P 2i (ik) , 

PD -i = t2 z2 ( 6 ) 

P 3 j (ik ) 2 + (ik) + P 4 j 

End j+E 21 & j etc. are the zero-lag response; and the variables a j and & j are 
defined as 

dj = ik (Xq-’ e^ kt 
& j = -k 2 a 0 j e ijkt ’ 

to be consistent with higher order terms. 

Finally, a general expression for the nonlinear aerodynamic models 
is constructed. This expression is to interpolate a complete set of harmonic- 
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oscillatory data with different frequencies to result in a formula for all 
frequencies. The response calculated at the test frequency encompasses the 
classical linear theory. It is noted that the response in time domain is given 
by the real part of the summation from each mode in Eq.(5). The 
expressions for all aerodynamic force and moment models are similar and 
will be determined with the following optimization procedures. 

2.2 Inner Optimization Procedure Without Constraints 

The main objective in the inner optimization procedure is to 
determine the E and H values for each mode independently by data-fitting 
the Aj and Bj at different frequencies. 

In the previous study, the first term in Eq.(5) was assumed to be 
constant based on a set of experimental data for a 70-deg. delta wdng. After 
two more sets of experimental data have been examined in this research, 
it is found that a linear polynomial function of k for A 0 (k) is better than a 
constant, especially in modeling C D responses. The Cj, which are given and 
unchanged in the inner procedure, are reference values and used to 
normalize the response given by Aj-iBj in the least-squared-error method. 

The inner optimization procedure mainly consists of two major 
methods. They are 

(l)Least-squared-error method: For fixed E and H values etc., the 
P 1 j,P 2 j,P 3 j,P 4 j are determined by minimizing the sum of squared errors. 
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(2)Gradient method: The values of E and H etc. are varied so that the 
sum of squared error is minimum. 

According to the previous study, a two-level gradient method is more 
effective than a straightforward gradient method because of the type of 
nonlinearity in this optimization problem. The detail about the least- 
squared-error method and gradient method are discussed in the following. 
The idea is illustrated in the flow chart of Fig.2. 

2.2.1 Least-Squared-Error Method 

At a specified reduced frequency, the magnitude of Fourier 
components in Eq.(4a) written in complex form is Aj-iBj. So 

A j- iB j = CjOC 0 j * [E x jik + E 2 j (-k 2 ) (7) 

+ (amplitude function) j* (1-PDj) ] 

For the assumed values of E and H etc., the only unknown variables in 
Eq.(7) are the coefficients of Pade approximants. Eq.(7) is rearranged as 
follows with all known values being on one side: 


i Aj - iBj - E xj ik - E 2 j <-k 2 ) 

(amplitude function) 

( 8 ) 

Pij (ik ) 2 + P 2j (ik) 

P 3 j (ik ) 2 + (ik) + P 4 j 
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It should be noted that Aj-iBj in Eq.(8) has been normalized by Cja 0 J . If both 
sides of Eq.(8) are multiplied by the denominator of the Pade approximant 
and separated into real and imaginary parts, Eq.(8) becomes 

Re = P^k 2 - PgjVjk 2 + P 4j Vj - Wj k = 0 (9a) 

and 

Im = P 2j k + P 3 jWjk 2 - P 4j Wj - Vjk = 0 (9b) 

Then, the coefficients of Pade approximants, i.e. P^ , are chosen by 
minimizing the sum of squared errors with different k’s. That is, 

Err = E Re(k i ) 2 + E Im (k ± ) 2 

By equating the first derivatives of Eq.(10) with respect to variables 
Pjj.Pgj.Pgj and P 4 j to zero, the coefficients of Pade approximants can be 
determined from the following equation: 

’ Ek i 
o 

-Evikj 
_EVik? 

where i is the index for different frequencies used in acquiring the 
experimental data, and the mode subscript j on V and W has been omitted. 


-Ev ± kJ 

Ew ± ki 


EVik 2 

-Ew i k i 


Ew ± kJ E(V 2 ki+W 2 kJ) -E(V?k?+W?k 2 ) p 3j 


EWiki E (v 2 k 2 +w 2 k 2 ) E(Vi+w^ 
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2.2.2 Gradient method without constraints 




Once the unknown coefficients P 1 j ) P 2 j,P 3 j and P 4 j are found, a one- 
dimensional gradient method is adopted to find better E and H values 
which minimize the sum of squared errors. The sum of squared errors is 
defined as : 

S = 52 (Aj) numerical 1 + 52 [Bj- (Bj) numerical! (12) 

where = Re[ RHS of Eq.(7) ] 

^numerical = Im [ RHS of E Q-(7) ] 

Because of the difficulty in locating a global minimum in a straightforward 
gradient method, a two-level gradient method described in the following is 
used in this investigation. 

(1) First level 

Step 1: The current E or H value is perturbed by a small amount AE 
or AH to find the gradient of sum of squared errors. 

Step 2: Each E or H value advances one step to new value according 
to its local gradient obtained in step 1. 

Step 3: Repeat Step 1 and Step 2 for each variable E or H until 
several iterations has been reached(it is set to be 5 in the current program) 

(2) Second level 

Repeat the first level until several iterations has been reached(it is 
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set to be 10 in the current program). 

2.3 Outer Optimization Procedure With Constraints 

As mentioned earlier in the introduction, the results for E and H 
values by the inner optimization procedure are not unique. They tend to 
vary with different initial guesses of E and H values. Moreover, not all of 
these results are acceptable due to the requirement that the denominators 
of Pad6 approximants must have real negative roots only, the latter 
requirement is necessary so that in time domain the corresponding 
exponential terms will die out at large time. A straightforward method 
which directly adds this requirement as constraints to the inner 
optimization loop has been tried. Unfortunately, for most cases, it was 
found that directly imposing constraints to the two-level gradient method 
tends to make the final results close to the initial guesses of E and H 
values. In other words, the constraints block the path of searching vector 
obtained in the gradient method so that a global minimum solution can not 
be reached in such a situation. To avoid this problem, an outer optimization 
loop has been added outside the inner optimization procedure. The outer 
optimization procedure is based on the same search technique as the inner 
one, i.e. the gradient method, except that additional constraints are 
imposed for the purpose of choosing proper initial values for E and H. A 
question arises as to what would be the objective function of the outer 
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optimization procedure. From the experience of validating the indicial 
formulation, an objective function which tends to minimize the lag effect of 
phase function in time domain at t=0 is chosen. That is, (l.-a lj" a 2p 1S 
minimized. 

A further examination of Eq.(5) is needed at this time. The remaining 
variables which are not calculated by any optimization procedure are Cj. 
Theoretically, E and H values are inversely proportional to Cj for a given 
response Aj-iBj, and can be adjusted to take appropriate values in the 
gradient method for a given reference value of Cj. But, numerical 
experimentation shows that this is not true. The reason is that the given 
response Aj-iBj in Eq.(8) is normalized by Cja 0 ^, while the E and H values 
always start from the same built-in initial guesses. In other words, the Cj 
must be properly chosen in a way to make the present method more 
workable and user-friendly in analyzing any given set of test data. 
Therefore, one more outer loop has been added outside the first outer 
gradient method to determine the best value of Cj within a feasible range 
of Cj. To illustrate the necessity of this additional outer optimization loop, 
test data for the 70-deg. delta wing(Ref.l3) used in Chin and Lan’s paper 
are reanalyzed by the modified method. Significant improvement at small 
time can be seen in the constant-rate pitch-up motion(Fig.3). The detail of 
the two-level outer optimization procedures are discussed in the following, 


M 
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and the flow chart is presented in Fig.4. 
2.3.1 Gradient Method With Constraints 




i m 


B 


B 


Unlike the gradient method in the inner optimization loop, the 
gradient method used in the outer loop is imposed with a constraint which 
requires all roots of the denominators of Pade approximants to be real and 
negative. In addition, another special constraint which needs to be applied 
only to the analysis of C D responses is that E n can not be negative 
according to the classical linear theory. By treating the whole inner 
optimization procedure as a "function box", a gradient method with 
constraints which comprises a gradient-search process and a series of 
checking-up has been implemented. This includes the following three steps. 

Step 1: Find an acceptable starting point 

For a given Cj, the starting initial guesses for E and H values are 
assumed to be a set of optimal initial values from the previous Cj or a set 
of build-in initial data in case the search of the previous Cj finds no 
reasonable starting values. If this starting point violates the constraints, 
each of E and H values will take turn to advance one step(maximal number 
of steps is set to be 4) until a reasonable starting point is found. 

Step 2: Construct an acceptable gradient vector 

Once an acceptable starting point has been found, a gradient vector 
is constructed by separately perturbing a small amount AE or AH for the 
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current initial guess. If the perturbation in any of E or H values leads to 
an unacceptable result, the gradient in the corresponding direction will be 
set to zero. It is noted that the objective function in this gradient method 
is defined as 

F^-l . / (1 . -a 1;j -a j2 ) (13) 

from numerical experiments, where a^ and a 2 j represent the lag effects at 
very small time. 

Step 3: Advance to a new acceptable point 

Based on the gradient vector in Step 2, a new set of initial values for 
E and H are obtained by advancing a given step length DS. Once again, the 
reasonableness of the new point should be checked. If it is an unacceptable 
point, the searching procedure must be stopped; and the last acceptable 
point found in Step 1 would be recalled to be the optimal set of initial 
values of E and H for the given Cj. Otherwise, Step 2 to Step 3 would be 
repeated further until several iterations has been reached(it is set to be 10 
iterations in the current program). 

2.3.2 Determination of reference value Cj 

In the search of the best reference value for Cj, no gradient method 
is needed. A simple way used in the present method is to evaluate different 
Cj over a feasible range. Based on the experience in analyzing the test data 



presented in this report and the earlier test data(Ref.l3) used by Chin and 
Lan, the optimal Cj usually satisfies the following inequality. 





Here, is in radian and j is the mode index. So, a set of feasible Cj for 
the second level search can be built up. Then, a cost function which is 
defined as 

F 2 = |a l3 + a 2j | + |a 1;j | + |a 2j | (14) 

is used in the evaluation of Cj. The best reference value of Cj is determined 
by selecting the Cj with a minimum cost function. 


m 
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3. 1NDICIAL FORMULATION 


3.1 General Mathematical Formulation 

According to Tobak(Ref. 14), in linear theory exponential functions 
which are obtained by applying the Fourier integral to the phase functions 
in frequency domain are good representations in time domain for 
aerodynamic response corresponding to a step change due to the asymptotic 
nature of exponential functions. In Chin and Lan’s paper this approach was 
applied to the nonlinear model and successfully set up a general 
mathematical formulation, called "indicial formulation". A brief review of 
this formulation will be given in the following. 

Once the appropriate coefficients in Eq.(5) have been found for a set 
of harmonic responses, Eq.(5) can be rewritten in a compact form by 
absorbing Cj into E and H values as follows: 

m 

C L = A 0 (k) + E i (E lj a j + E 2 j ftj) (15) 

m 

+ E (amplitude function) j * (phase function) j 
j=i 

where phase functions are defined as 
1 - p ij < ik > 2 + P 2 j < ik > 

p 3j (Ik)* + Ik + p 4ji (16) 

q _ i(jk) a x j _ i(Dk) a 2j Ufc>; 

i ( jk) + ja 3 j i ( jk) + ja 4 j 
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Based on a step input (Ref.4 ), the phase function will be converted into 
exponential functions in time domain but the amplitude function is kept 
unchanged. Therefore, an indicial response C L for the circulatory lift, i.e. 
the third term in Eq.(15), can be defined as 


C 


^indicial 


m 

£ (amplitude function) 


(l.-a lje " a ^ jt -a 2 j e- a ^ jt ) 


(17) 


From Eq.(l), it is known that aerodynamic response at time t is made of 
three parts. The first part is the initial response at the starting point. The 
second part is the response without hysteresis effect; while the last part is 
the indicial response due to changes in a and a in the prior motion. The 
initial response C L (0) is obtained by setting the running variable x to zero 
in the indicial response C L and the zero-lag response term is identical in 
every harmonic motion. Then, the integrands in the third part can be 
determined explicitly by taking derivatives with respect to a and a from 
Eq.(17). It should be noted that an average value C ave also will be added in 
Eq.(l) because a in the amplitude functions of Eq.(17) denotes only a 
perturbation from o^ in the harmonic model. Therefore, the final form of 
indicial formulation for arbitrary motions is given as 
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C-L(t') - C Lindiclal 


[t / -T , OC(X ) , OC (T ) ] x=0 + c 

ave + 


m c t ! d (a. f . ) 


■* ( 1 -a x je 


m /*, 

E f 1 

j-iJo 

— E f t7 lii±il*(i-a 1 . 
v^iJo do 


, ft _ a 3 j j (t / -x) ^ e -a 4 -j j ( t / —x ) j da(x) 


dx 


dx 


-a 3;) j (t'-T) _ a ^ e -a 4;) j (t'-T) ) da(X ) 


dx 


dx 


(18; 


3.2 Implementations of Indicial Formulation 

To perform the time integration in Eq.(18) for arbitrary motions, 
some special numerical implementations have been made in the present 
program. First of all, a 3-point Simpson rule is chosen to evaluate 
numerical integration. Secondly, an arbitrary motion has to be represented 
locally by a cosine function in order to utilize the results of harmonic 
modeling. At a certain time of an arbitrary motion, the total angle of attack 
o^ and ct 1 can be described by the cosine and sine functions as 

cq (x ) =a m +a Q cos (kx +<J> ) (19a) 

Oq ( X ) = - 0t o ksin(kx+(j)) (19b) 


20 


By knowing the harmonic model’s mean angle of attack ot,^ and amplitude 
a 0 , an equivalent frequency k and an equivalent phase angle $ at a given 
instantaneous time x can be solved by Newton’s method. To smooth out 
possible discontinuity in response when the given motion has a sudden 
change in a , an a 0 which is slightly greater than the actual amplitude is 
frequently used. It should be emphasized that this dose not change the 
instantaneous values of a x and d-L in the actual time history of motion. It 
merely changes the values of equivalent frequency k and phase function <)). 
Finally, the C ave in Eq.(18) which was assumed to be a constant in the 
previous study has been reformulated. A constant term Aq in the harmonic 
modeling is now replaced with a linear polynomial function of k for Ao(k). 
Based on the similar concept in the indicial response, the C ave at a given 
t’ should be made of all Ao(k) of prior states in the motion up to the current 
time. So, a simple way to calculate C ave is to take an average, i.e. 

Cave = 1 fiA 0 (k m )]/M ( 20 ) 

m=l 

where m is an index for running time, and k m is the equivalent frequency 
based on the original amplitude in the harmonic model. 

In addition to the above-mentioned general implementations in the 
indicial formulation, several special treatments used in some particular 
motions will be discussed in the next chapter. 
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4. RESULTS AND DISCUSSION 
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4.1 Test Data 

In the present study, test data for three kinds of motions obtained at 
the NASA Langley Research Center are taken to be analyzed or compared 
with numerical results calculated from simulations of the specified motion. 
They are described in detail as follows: 

(1) Large Amplitude Harmonic Motions: Test data of harmonic 
motions with large amplitude are used for numerical modeling in Fourier 
functional analysis. In this category, there are two groups of test data, and 
each group includes three sets of aerodynamic data for 0^,0^ and Cj^ with 
five reduced frequencies. The first group is for a 70-deg. delta wing, and the 
second one is for an F-18 model. Both test models are harmonically 
oscillated about 32.5 deg. of angle of attack with an amplitude of 30 degree. 
This large coverage in angle of attack which contains low a region, near- 
stall region and post-stall region exhibits a good characterization of 
hysteresis behaviors. In addition, the test data for a 70-deg. delta 
wing(Ref.l3) used in Chin and Lan’s paper also have been reanalyzed for 
the purpose of comparison. 

(2) Constant-Rate Pitching Motions: Test data of constant-rate 
pitching motions are used for validation of the indicial formulation in 
arbitrary motions. The experimental data in this category were taken from 
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the same models as described in the large amplitude harmonic motions. 
Responses with three pitch rates for each model have been measured in 
pitch-up motions as well as pitch-down motions. The angle of attack starts 
from 2.5 degree and stops at 62.5 degree in the pitch-up motions; but moves 
in reverse in the pitch-down motions. 

(3)Medium Amplitude Harmonic Motions: An additional group of 
forced harmonic oscillation data for an F-18 model provides a further 
examination of the interpolative property of the mathematical model. The 
mean angle of attack and amplitude of oscillation in this harmonic motion 
are intentionally set to be different from those in the large amplitude 
harmonic motions. In the present test, the F-18 model was harmonically 
oscillated about 22.5 degree angle-of-attack with an amplitude of 20 degree. 
4.2 Fourier Functional Analysis 

As pointed out in chapter 2, a constant term for A 0 (k) in the previous 
study is replaced by a linear polynomial function of k, especially while 
analyzing C D responses. To show the necessity of this change, the term A 0 
versus the test reduced frequencies for the F-18 C D responses are plotted 
in Fig. 5(a), and also the numerical results by both models are compared 
with experimental data at the highest k in Fig.5(b). Both figures show that 
a discrepancy of about 0.1 exists if an average value for A^Ck) is applied. 
Regardless of the optimization procedures, other task to improve the 
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accuracy of the present method is to include the static test data as 
additional dynamic stall data at a very low reduced frequency, such as 
k=1.0E-6. The main purpose of this implementation is to avoid possible poor 
extrapolation at low reduced frequency. 

The results by the Fourier functional analysis are presented in Fig.6 
for the 70-deg. delta wing and Fig.7 for the F-18 model respectively. Five 
Fourier terms are used in this analysis and the calculated coefficients for 
the two aerodynamic models are listed in Table(l) and Table(2) respectively. 
All results were done by the same set of build-in initial data for E, H and 
Cj without any try-and-error guess by a user. The results show that the 
present method is able to capture all major hysteresis effects. Compared 
with the up-stroke data, most of the down-stroke data are modeled with 
less accuracy. The reason can be found by taking a further look at the C L 
test data for the presently-used 70-deg. delta wing (see Fig. 8(a)) and those 
for the other 70-deg. delta wing(Ref.l3) used in Chin and Lan’s paper (see 
Fig.8(b)). In Fig.8(b), it appears that the flow becomes attached in the low 
a region (0-20 deg.) during the down-strokes. On the contrary, strong 
hysteresis effects still exist under similar conditions in Fig. 8(a). So, the 
trend of the hysteresis behavior on down-strokes is not as consistent as that 
on up-strokes in the presently-used test data. This may imply that a higher 
order Pade approximant could be needed to model responses which have 

24 


more complicated hysteresis effects. It is noted that the mismatched part 
of F-18 C M responses at k=0.015(see Fig.7(c)) is due to the even more 
inconsistent trend of the hysteresis behavior from one frequency to the 
other occurring in the region of high angle of attack. 

4.3 Indicial Formulation 

4.3.1 Harmonic Motion and Harmonic Ramp Motion 

The harmonic motion is used to compare with Fourier modeling 
results which have been well fitted with test data. On the other hand, the 
harmonic ramp motion can be used to show how good the agreement with 
the static value is at the time when the motion stops. 

As mentioned in chapter 2, discontinuity could happen in the 
calculation of time integral if the given motion has a sudden change in a . 
This can be easily solved by slightly increasing the amplitude of the original 
model, such as by 2.5 degree. In the following calculations, the amplitude 
(ocq) which was 30 degree in the testing is set to be 30.5 degree for the 
harmonic ramp calculation and 32.5 degree for the constant-rate pitching 
motions respectively. 

The results of harmonic motions and harmonic ramp motions are 
plotted in Fig.9 for the delta wing and Fig.10 for the F-18 model 
respectively. In the harmonic ramp motions, all responses eventually 
approach the static value corresponding to the angle of attack when the 
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motion stops. 

^ j ;; 

^ 4.3.2 Constant-Rate Pitching Motion 

- • This is used to compare with test data at the same pitch-rate. 

Two special treatments about this type of motion should be 

— mentioned here. First of all, a question arises about the results by Eq.(19), 
when the actual angle of attack (a 1 in Eq.(19)) is near the two ends of a 
harmonic model’s a range, for example 2.5 or 62.5 deg. in the present test 

| 1 ........ 7 ■ 

— model. In such a situation, the equivalent k obtained from Eq.(19) tends to 
be high because the corresponding d x is too large comparing with a in the 

V 

harmonic model. From the experience in calculating the constant-rate 
“ pitching motions, it was found that an unreasonably extrapolated high 

1 value of k at a starting point would lead to an unacceptable result in 

simulation. So, one of the variables, a m and a 0 , must be treated as an 
unknown instead of k, when the extrapolated k-value is greater than a 
^ given allowable value k max . Through a series of tests on both variables, the 

amplitude of harmonic model a 0 was chosen as the other unknown in case 
the extrapolated k-value exceeds the given allowable k max . Therefore, 
_ ' Eq.(19) will be replaced by the following equation if the equivalent 

frequency k in Eq.(19) is larger than the given k max : 
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a i (x ) =a m +a 0 cos (k max T +<])) 


(21a) 


Cci (x) =-a 0 k max sin (k max x+(t)) (21b) 

where a 0 and <]) are the unknowns, and are, again, solved by Newton’s 
method. 

Secondly, as the constant-rate pitch-down motions start at a high 
angle of attack, the time integration should start from a static value by 
setting t’— to the first term of Eq.(18). 

The results of the constant-rate pitching motions are presented in 
Fig. 11 for the delta wing and Fig, 12 for the F-18 model respectively. As 
expected, all results for pitch-up motions are well predicted except in the 
region near the starting point. The reason for this is that the phase lag 
terms (l.-a 1 j-a 2 j) for x=0 in the first term of Eq.(17) are not able to perfectly 
represent a starting point which physically does not involve any hysteresis 
effects. On the other hand, some of the results for pitch-down motions are 
not as good as those for pitch-up motions. There are two possible reasons 
for this. Firstly, the harmonic model which does contain large hysteresis 
effects in the high angle-of-attack region can not precisely depict an initial 
response at the starting point in the high a region which definitely 
possesses a static value without any phase lag. Secondly, the poor 
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numerical modeling of harmonic motions in down-strokes which is caused 
primarily by the non-smooth variation of the response from one frequency 
to the other, especially at low and moderate k, also is a reason for these 
discrepancies. 

4.3.3 Validation of the Interpolative Property 

Harmonic responses with a lower mean angle of attack(a m ) and a 
smaller amplitude(a 0 ) for the F-18 model are calculated by using the 
indicial formulation to compare with the test data in the same conditions. 
The aerodynamic models are still those based on the test data with 0^=32. 5 
degree and ocq= 30 degree. 

Since the a range and reduced frequencies still are within the range 
of the original harmonic model, no extrapolated equivalent reduced 
frequencies are expected. For the known otj and ctj time histories, the 
equivalent k and 4> at certain time % can obtained from Eq.(19) without any 
change in a 0 and of the original harmonic model. To avoid the error 
which usually happens at a small time, the integral in Eq.(18) are 
evaluated from the third cycle and only over one period for the periodical 
motion. 

The results are plotted in Fig. 13. Comparing with the experimental 
data, the simulation by using the indicial formulation shows reasonably 
good agreement. It is also shown that the higher the reduced frequency is, 
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the better the calculated responses are. The reason is that the simulation 
for higher reduced frequency includes higher equivalent k which are well 
modeled in the original harmonic modeling than the simulation for lower 
reduced frequency. Moreover, some other results should also be noted in 
this simulation, first of all, the results of Cp responses reveal that using a 
model containing large hysteresis effects to calculate a motion with less 
dynamic effect may not be as good as expectedfsee Fig. 13(b)). Secondly, 
although the analysis for C M responses fails to model accurately the high 
a region (a>40 deg.) at low and moderate reduced frequencies (see Fig.7(c)), 
the harmonic motions with a oc-range below 42 deg. still can be well 
simulated by the indicial formulation. 


5. CONCLUSIONS 


In the present research, a method involving Fourier functional 
analysis and indicial formulation proposed by Chin and Lan has been 
extensively examined and modified for the purpose of applications to 
airplane configurations. Two sets of test data from NASA Langley Research 
Center which include a 70-deg. delta wing and an F-18 model have been 
used to show the applicability of Fourier functional analysis and validate 
the indicial formulation. Extensive examination of this method has led to 
the modification of the method by adding an outer optimization procedure 
with constraints to the original optimization loop to automatically choose 
appropriate starting values for the unknowns in the nonlinear optimization 
process. 

The results from the Fourier functional analysis showed that the 
general expression for the aerodynamic response to harmonic motions 
throughout a range of k was capable of accurately modeling nonlinear 
responses with large phase lag except in the region with an inconsistent 
hysteresis behavior from one frequency to the other, such as the F-18 C M 
responses at low reduced frequencies. The results by time integration of the 
indicial integral showed the applicability of the latter to harmonic motions 
and ramp-type motions. Moreover, all pitch-up motions with constant rate 
were well simulated by indicial formulation; while the results for the 
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corresponding pitch-down motions were produced with less accuracy but 
still in a correct trend. Finally, The results for the F-18 model’s 
aerodynamic response to harmonic motions with different mean angle of 
attack and amplitude was indicative of good interpolative property of the 
present aerodynamic models. 
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(a) C L Responses 
A 0 (k)=Q. 567+0. 140*k 
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H 5j 
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0.850 


II 
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0.377 
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■H 
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EH 
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(b) C D Responses 
A o (k)=0.476+0.8248*k 
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2 
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-11.78 
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-.0538 
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4 

12.31 
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18.868 
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5 
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1.084 

-0.677 

-.0010 
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| 


Table(l) Model Coefficients for A 70-deg. Delta Wing. 
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(c) C M Responses 
A 0 (k)=0.046+0.Q773*k 
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Table(l) Concluded. 
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(a) C L Responses 
A 0 (k) = 1.131+0.804*k 
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(b) C D Responses 
A 0 (k)=0.916+1.552*k 
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Table(2) Model Coefficients for An F-18 Model. 
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(c) C M Responses 
A 0 (k)=-0. 07614+0. 7105 *k 


j 

C i 

EH 


H H 

H 2j 

H 3j 

H 4i 

H 5] 

EH 

1 

2.5 

-0.942 

-0.954 

-0.150 

-0.850 





2 

BE 

-0.098 

-0.269 

-0.350 

0.780 

-1.001 




3 

1.75 

-0.117 

-0.112 

-0.250 

-0.050 

0.250 

-0.050 



4 

3.0 

-0.098 

0.021 

-0.050 

0.798 

0.994 

-0.056 

-0.050 


5 

23.0 

-0.018 

-0.005 

0.050 

0.856 

-0.066 

-0.050 

-0.050 

-0.05 

j 

P H 

EH 

EH 

ESI 

a ii 

a 2i 

a 3i 

.. a 4j 


1 

-1.464 

0.667 

1.015 

0.001 

0.6702 

-2.112 

-.0010 

-.9840 


2 

22.400 

0.141 

23.931 

0.001 

0.1241 

0.8119 

-.0010 

-.0408 


3 

35.605 

0.277 

26.874 

0.001 

0.2543 

1.0705 

-.0010 

-.0362 


4 

4.007 

1.086 

5.440 

0.018 

1.2947 

-0.558 

-.0207 

-.1631 


5 

3.289 

0.926 

4.790 

0.001 

0.9316 

-0.245 

-.0010 

-.2077 



Table(2) Concluded. 
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Figure 1 A New Air Combat Tactic with Post-Stall 
Maneuvering 


38 




Figure 2 Flow Chart for The Inner Optimization 
Procedure. 
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Original Method 
Modified Method 


Comparison of Unsteady Lift Coefficient in 
Constant-Rate Pitch-up Motion by The Modified 
Method with The Results by The Original 
Method for A 70-deg. Delta Wing. 





<define> 1 0 P = inner optimization procedure (fig. 2) 



Figure 4 Flow Chart for The Outer Optimization 
Procedure. 
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Figure 6 Continued. 
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Figure 6 Continued. 
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Figure 6 Concluded. 
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Figure 7 Continued. 
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Figure 10 Concluded. 
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Figure 11 Continued. 







(b)Drag Data 


i — i — i — i — j — i — i — i — i — | — i — i — i — i — | — i — i — r 

pitch-down 



i i L_ i _i i i l 1 l l i _i 1 1 I 1 

) 30 40 50 

a 

(b)Drag Data 
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Figure 12 Continued. 
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Figure 12 Concluded. 
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